/**********************************************************
read me:
	1. n can be up to 33,553,920(testing)
	2. haven't implemented scan algorithm
**********************************************************/

#include <stdio.h>
#include <sys/resource.h>

#include "body.h"
#include "cutil_math.h"

#define spaceDim 3

#define BLOCK_WIDTH 512

#define G 6.67384e-11f

#define MAX_NUM_STARS 33553920

__constant__ float softeningSquared =0; //the square of the softenting factor->EPS

/****************************************************************************************************
for calculate accelerations
****************************************************************************************************/
//__device__ body-body
__device__ float3 
bodyBodyInteraction(float3 ai, float4 bi, float4 bj) 
{
    float3 r;

    // r_ij  [3 FLOPS]
    r = make_float3(bj)-make_float3(bi);

    // distSqr = dot(r_ij, r_ij) + EPS^2  [6 FLOPS]
    float distSqr = r.x * r.x + r.y * r.y + r.z * r.z;
    distSqr += softeningSquared;

    float distSixth = distSqr * distSqr * distSqr;
    
    //for current case only-> soften factor is 0, when it is not 0 then don't need this if
    if(distSixth!=0){
	    float invDistCube = 1.0f / sqrtf(distSixth);
	    
	    ai += (G*bj.w*invDistCube)*r;;
    
    }
    

    return ai;
}

//__device__ tile
__device__ float3
tile_calculation(float4 myPosition, float3 accel, float4* shPosition, unsigned int tileId, bool rest, unsigned int restNum)
{
	int i =0;
	//last not copleted filled block/tile 
	if(rest&&tileId==(gridDim.x-1)){
		for(i = 0; i < restNum; i++){
			accel = bodyBodyInteraction(accel, myPosition, shPosition[i]);
		}
	}
	//fully filled block/tile
	else{
		for(i = 0; i < blockDim.x; i++){
			accel = bodyBodyInteraction(accel, myPosition, shPosition[i]);
		}
	}
	
	return accel;
}

//__global__ do the calculation
__global__ void 
updateVel_shMem(void* devX, void* devA, int N, bool rest, unsigned int restNum)
{
	__shared__ float4 shPosition[BLOCK_WIDTH];
	
	float4* globalX = (float4* )devX;
	float4* globalA = (float4*)devA;
	float4 myPosition;
	int i =0;
	int tile = 0;
	float3 acc;// = {0.0f, 0.0f, 0.0f};
	int gtid = blockIdx.x*blockDim.x + threadIdx.x;
	
	
	myPosition = globalX[gtid];

	acc = make_float3(globalA[gtid]);

	for(i = 0, tile = 0; i < N; i +=blockDim.x, tile++){
		
		if(gtid < N){
			int idx = tile*blockDim.x + threadIdx.x;
			if(idx < N){
				shPosition[threadIdx.x] = globalX[idx];
			}
		}
		
		
		__syncthreads();
	 	
	 	if(gtid < N){
	 		acc = tile_calculation(myPosition, acc, shPosition, i, rest, restNum);
	 	}
		
		__syncthreads();
	}
	//save the result in global memory for the integration step
	if(gtid < N){
		float4 acc4 = make_float4(acc);
		globalA[gtid] = acc4;
	}
	
	
	
}

/****************************************************************************************************
update positions
*****************************************************************************************************/
__global__ void updatePos(float4* pos, float4* vel, int n, float dt)
{
	int x = blockIdx.x*blockDim.x + threadIdx.x;

	if (x < n)
		pos[x] += dt*vel[x]; // 3 flops
}
/****************************************************************************************************
test: print of time
*****************************************************************************************************/
void
print_times(const char *kernel, struct rusage *before, struct rusage *after) {
  printf("---- %s -----\n", kernel);
  printf("[sys a] sec: %ld usec: %ld\n", after->ru_stime.tv_sec, after->ru_stime.tv_usec);
  printf("[sys b] sec: %ld usec: %ld\n", before->ru_stime.tv_sec, before->ru_stime.tv_usec);
  printf("[usr a] sec: %ld usec: %ld\n", after->ru_utime.tv_sec, after->ru_utime.tv_usec);
  printf("[usr b] sec: %ld usec: %ld\n", before->ru_utime.tv_sec, before->ru_utime.tv_usec);
  printf("sys: %ld usr: %ld\n", after->ru_stime.tv_sec - before->ru_stime.tv_sec, after->ru_utime.tv_sec - before->ru_utime.tv_sec);
  printf("\n");
  fflush(stdout);
}

/****************************************************************************************************
main
*****************************************************************************************************/
int
main(void){
	//test time
	struct rusage times[2];
	
	//cpu
	bodyList* bl = readBodyList();
	unsigned int n = bl->size;
	bool rest = false;
	unsigned int restNum = 0;
	
	//gpu
	float4* d_pos =NULL;
	float4* d_vel = NULL;
	
	//allocate memories on gpu
	d_pos = safe_cudaMalloc(n, bl);
	d_vel = safe_cudaMalloc(n, bl);
	
	//memcpy to device
   	size_t size = n*sizeof(float4);
   	cudaMemcpy(d_pos, bl->pos, size, cudaMemcpyHostToDevice);
	cudaMemcpy(d_vel, bl->vel, size, cudaMemcpyHostToDevice);
   	
   	//calculate grid width and see if there needs an extra block
   	unsigned int GRID_WIDTH = n/BLOCK_WIDTH;
   	restNum = n%BLOCK_WIDTH;
   	if(restNum!=0){
   		GRID_WIDTH += 1;
   		rest = true;
   	}
   	
   	//determine the grid and block size
   	dim3 dimb(BLOCK_WIDTH, 1);
	dim3 dimg(GRID_WIDTH, 1);
	
	//invoke kernel2 with shared memory
	getrusage(RUSAGE_SELF, &times[0]);//timestamp
	
   	updateVel_shMem<<<dimg, dimb>>>((void*)d_pos, (void*)d_vel, n, rest, restNum);
   	cudaThreadSynchronize();
   	
   	getrusage(RUSAGE_SELF, &times[1]);//timestamp
	print_times("shared_memory", &times[0], &times[1]);
   	
   	//invoke kernel to update the positions
	updatePos<<<dimg, dimb>>>(d_pos, d_vel, bl->size, 1.0f);
	
	
   	
   	//copy accels back to cpu
   	cudaMemcpy(bl->pos, d_pos, size, cudaMemcpyDeviceToHost);
	cudaMemcpy(bl->vel, d_vel, size, cudaMemcpyDeviceToHost);
   	
   	// print output
	/*int i;
	for (i = 0; i != bl->size; ++i)
	{
		float4* pos = &bl->pos[i];
		float4* vel = &bl->vel[i];
		printf("%9.2f %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f\n", pos->x, pos->y, pos->z,
		//printf("%9.2e %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e\n", pos->x, pos->y, pos->z,
		vel->x, vel->y, vel->z, pos->w);
	}*/
   	
	//free gpu
	cudaFree(d_pos);
	cudaFree(d_vel);
	//free cpu
	freeBodyList(bl);
	
	return 0;
}
